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Abstract: The Phase-Locked Loop (PLL) is a key component of modern electronic 
communication and control systems. PLL is designed to extract signals from transmission 
channels. It plays an important role in systems where it is required to estimate the phase 
of a received signal, such as carrier tracking from global positioning system satellites. In 
order to robustly provide centimeter-level accuracy, it is crucial for the PLL to estimate the 
instantaneous phase of an incoming signal which is usually buried in random noise or some 
type of interference. This paper presents an approach that utilizes the recent development 
in the semi-definite programming and sum-of-squares field. A Lyapunov function will be 
searched as the certificate of the pull-in range of the PLL system. Moreover, a polynomial 
design procedure is proposed to further refine the controller parameters for system response 
away from the equilibrium point. Several simulation results as well as an experiment result 
are provided to show the effectiveness of this approach. 

Keywords: non-linear systems; phase-locked loop; optimization 
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1. Introduction 

For nonlinear control systems, one would often like to know the domain of attraction (DoA) of an 
equilibrium point. Often, this domain is difficult to both find and represent computationally. The usual 
mathematical tool used for analyzing of the region of attraction is Lyapunov's method. This gives us 
a sufficient condition for local stability, although it is often difficult to find a Lyapunov function that 
can be used as a certificate for the whole DoA. Several prior approaches have used quadratic functions, 
for example [1—3]. In particular, the approach of [3] makes use of semidefinite programming to find 
a quadratic function whose sublevel-set is a good inner approximation to the region of attraction. For 
system in which the DoA is complicated, an ellipsoid may not provide a good approximation, and the 
above methods leave a large unexplored region within the DoA. 

With recent developments in algebra and sum-of- squares techniques, it is now possible to solve for 
a Lyapunov function with a more general polynomial form [4,5]. Positive definiteness properties are 
replaced by sum-of-squares(SOS) constraints which can be efficiently solved using convex optimization. 
The SOSTOOLS [6] toolbox for MATLAB has been developed as an easy computational tool to solve 
problems that utilizes the SOS techniques. This approach has also allowed finding a Lyapunov function 
within some specified semi- algebraic region [7,8]. While Lyapunov approach provides a method to 
certify a given inner approximation to the domain of attraction, it does not immediately provide a way to 
find it. Tan [9] later extended this concept by using unions of SOS polynomials to estimate the domain of 
attraction but Tan's approach needs to solve a bilinear optimization. The level-set method [10] has been 
developed to find a semi-algebraic representation of the DoA by semidefinite programming. With these 
polynomial techniques, it is possible to precisely estimate the DoA of a nonlinear polynomial system 
and to find a suitable Lyapunov function as the stability certificate. 

A phase-locked loop (PLL) system is a nonlinear system with limited domain of attraction. Due to 
its importance in communication systems, analyzing and designing a PLL system has attracted many 
attentions in this field [11-17]. The current approach for designing a controller for a PLL system is 
still largely based on the linear model [17]. Hence, the performance of the resultant system cannot be 
guaranteed at system states far away from the designed equilibrium point. 

In this paper, we utilize the current SOS techniques to analyze the domain-of-attraction of a PLL 
system. A local Lyapunov function can then be found as the certificate of the DoA using the proposed 
approach. The Lyapunov function will be further used to improve the stability region and performance of 
the PLL system. Using this approach, it is possible to design a PLL with a predefined form of controller 
that has larger domain-of-attraction than the linear design approach. Moreover, since we are designing 
the system in the nonlinear region, system dynamics outside the linear region can be further refined. 
Examples of second order PLL systems are used later in this paper to show the effectiveness of this 
design approach. In the provided examples, we demonstrated that the system designed by the proposed 
method has 20% larger domain-of-attraction and less overshoot with faster convergent rate than the linear 
design approach. 

This paper is organized as follows. Section 2 contains some preliminary knowledge about SOS 
techniques. The advection algorithms as well as the SOS methods for finding a local Lyapunov function 
are stated in Section 3. A brief discussion of the configuration of a PLL can be found in Section 4 and 
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the proposed method for analyzing and designing a PLL controller are presented in Section 5. Then the 
paper is summarized in Section 6. 

2. Preliminaries 

The following are some definitions that will be used frequently in this paper. M[x] is used to 
represent the ring of polynomials in x with real coefficients. A polynomial / E M[x] is called positive 
semidefinite (PSD) if f(x) > 0, for all x E W 1 . A polynomial / is called SOS if there exist polynomials 

g h g s E R[x] such that / = gl + g%-\ Vg\. Clearly if / is SOS then / is PSD. It is also well-known 

that the converse is not true. £ denotes the set of all SOS polynomials in M[x], M + is used to represent 
the set of nonnegative real numbers. B r is used to represent the open ball with radius r centered at 
the origin. 

Suppose g : MJ 1 —> K is C l . Define the 0-sub-level set of g to be sub(g) C W 1 given by 
sub(g) = {16 1™ g{x) < 0 }. Further define the boundary of sub(g) as d sub(g). 

One feature of the proposed advection algorithm is that the advection problem can be converted into 
a semidefinite program. The following is a standard form of a semidefinite program. 

min trace (CX) 
x 

s. t. trace(AjX) = bi for i = 1, . . . , m 

XhO, 

where X E M. nxn is symmetric. X y 0 means that z T Xz is positive semidefinite for all z E W 1 . 

The condition of one semi-algebraic set containing another semi-algebraic set is one of the key 
constraints used in this paper. The following lemma shows that this kind of relationship can be converted 
to constraints on the coefficients of the polynomials. The proof can be found in [4] or [7]. 

Lemma 1. Given p, q E M[x], suppose there exist s 0 ? si 6 E such that 

s 0 -s iq + p = 0 forallxEW 1 (1) 

Then sub(q) C sub{p). Further, given q and the degree bound ofp, Sq, and s\, the set of coefficients of 
p, s 0 and Si satisfying (I) is the feasible set of a semidefinite program. 

Proof. See, for example, [4] or [7]. 

The representation shown in Lemma 1 is one of the simplest cases of Schmudgen's Theorem [18]. 
Schmiidgen's Theorem states that if p E M[x] is strictly positive inside a compact semi-algebraic set S 
generated by pi, . . . ,p m as S — {x E M. n \ pi > 0, i — 1, 2, . . . , m}, then 

p = ^ v p v i . ..p v ™s v 

where v = (vi, . . . ,v m ) E {0, l} m and s v E S. Putinar [19] later showed that under some additional 
constraints on pi, p has a simpler representation as 



V = s 0 + sipi + . . . + s m p m 
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The gap between Schmiidgen's and Putinar's representation is later investigated by Jacobi and 
Prestel [20]. In the simple case shown in Lemma 1, if sub(q) C sub(p) and sub(q) is compact, the 
representation of p by (1) is always possible. 

The following result is similar. Given q E M[x], if there exists s 0 , Si E S and e > 0 such that 

So + — P + e = 0 

then sub(p) C sub(q). 

Usually q is a given polynomial and p is the solution to find such that sub(p) and sub(q) approximately 
represent the same set with some other constraints on p, such as having lower degree or passing through 
several pre-specified points. The above results are used to construct such constraints. 

3. Acquiring the Local Lyapunov Function 

Finding a local Lyapunov function is coupled with finding the DoA. Without a clear knowledge of the 
actual shape of the DoA, it is hard to find a Lyapunov function that can be used to represent the entire 
DoA [4,8]. To deal with this difficulty, we utilize the current development in set advection [10]. 

3.1. Set Advection 



In this paper, we will consider the following autonomous system 

x(t) = f(x) (2) 

where / : R n — > R n is locally Lipschitz. From the basic local existence and uniqueness theorem [21], 
given an open subset U E W 1 , there exist c G R+ such that the autonomous system (2) has a unique 
solution for any z G U in the compact time interval [— c, c] . 

We define the flow map <f> t : R n x R — > R n to be the local unique solution of 

= /(&(*)) for t G [-c,c], c(z) GR +1 zE W 1 

0oO) = Z 

For any t G R such that <j>t{x) exists, the map 4> t '■ R n — > R n is continuous, invertible and has a 
continuous inverse, i.e., it is a topological homeomorphism on MJ 1 [22]. 

Given t E R, we define the time t advection operator A t : C(R n , R) — )• C(R n , R) by 

q = A t p if q(x) = p (0_ t (x)) for all i6R" 

where C(X, Y) is the set of functions mapping from X to Y. The map A t is also called the Liouville 
operator associated with /. A very important property is that it is linear. Figure 1 shows the concept of the 
advection operator. Given polynomial p, A t maps the coefficients of p to another polynomial q such that 
sub(q) = 4> t sub(p). We relate the advection operator to the advection of sets in the following remark. 

Remark 1. Suppose gi, g<i are functions mapping R n to R. If ' g% = A t gi then sub{g2) = (f>t (sub(gi)). 
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Figure 1. The advection operator A t . 




3.2. Time-Stepping 

Since we are performing advection, we must use an approximation to the flow map <p h with time step 
h. Many such approximations are possible, and are provided by the theory of numerical integration. The 
first-order Taylor approximation to q = A h p is the map B h : C(R. n , M) — > C(M. n , K) given by 

q = BhP if q(x) = p(x) — hDp(x) f (x) 

where the derivative Dp(x) is a linear map Dp(x) : W l — > MJ 1 at each point x. 

Based on the required accuracy of the advection, we could also choose to use higher order Taylor 
approximation. However, depending on the system dynamics, this usually will lead to the requirement 
of using higher degree polynomials in the sum-of-squares constraints. The relationship between the 
accuracy and the degree of polynomials will be further investigated in future work. 

3.3. Domain-of -Attraction Estimation 

The set advection concept is used to estimate the DoA of a system. We use the following definition 
of the DoA in this paper. 

Definition 1. Suppose f : W 1 — > M. n is analytic with the flow map, 0, and the origin is asymptotically 
stable. Define the domain-of-attraction (also called the basin/region of attraction) of f to be R C M. n 
such that for any x G R, <pt{x) is defined for allt > 0 and lim^oo <pt{x) = 0. 

The following properties can be easily derived. The detailed proofs can be found in [10]. 

Lemma 2. Suppose f is analytic and the origin is asymptotically stable and R ^ 0. Suppose also 
Si C R and 0 G Si, and Si is a connected closed positively invariant set. Let h > 0 be a positive 
constant, and define the backward advection of Si to be S 2 , given by 

S 2 = 4>-hSi 

Then Si C S 2 C R, and S 2 is also connected, closed and positively invariant. Furthermore, 

dS 2 = <j>- h dS h 
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Theorem 1. Suppose f is analytic and the origin is asymptotically stable and h > 0. Also suppose 
0 E So and So C R is a closed connected positively invariant set, such that there exists e > 0 such that 

Be C S 0 . 

Define the sequence of sets S 0 , Si, S 2 , ■ ■ ■ by 

S k+ i = (p^ h S k for k = 0,1,2, ... 
Then this sequence converges to R in the following sense: 

(i) S k C Rfor all keK 

(ii) Sk C Sk+ifor all k EN. 

( Hi) For every x E R, there exists n such that x E S n 

3.4. Star-Shaped Constraint 

For the case of estimating the DoA, we introduce the concept of star-shaped sets. The star-shaped sets 
have many important properties and can be easily implemented as a semidefinite program. We now start 
with the first property. The detailed information about the star- shaped set can be found in [10]. 

Definition 2. A set S E W 1 is called star-shaped if for all x E S 

\x ES for all A E [0, 1] 
The set S is called strictly star-shaped if for all x E S 

Xx E int(S) for all A E [0, 1) 

Note that a star shaped set S is connected. We now give a simple sufficient condition that ensures a 
sub-level set is star-shaped. We make the following definition. 

Definition 3. Suppose g : M. n — > R. We call g strictly star-shaped if g is C 1 and further satisfies 
g(0) < 0 and 

Dg(x)x > 0 for all x ^ 0 

The following lemma shows the connection between strictly star-shaped functions and 
star- shaped sets. 

Lemma 3. Suppose (j:l n 4lis strictly star-shaped. Then sub(g) is strictly star-shaped. 
Proof. Suppose x E sub(g). Let y : R + — > M. n be the function 

y (t) = e~*x 

The trajectory of y{t) follows the straight line connecting x and the origin. We would like to show that 
y(t) E sub(g) for all t > 0. We have 

j t g(y(t)) = -Dg(y(t))y(t) 
< 0 
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for all t > 0. Also 




and since y(0) = x we have g (y(t)) < 0 for alH > 0 as desired. 

Now we need to show that iiy E d sub(g) then there does not exist A E [0, 1) such that \y E d sub(g). 
Suppose for the sake of a contradiction that there does exist such a y and A. We know that there exists 
such A > 0, since g(0) < 0. Define the function h : [0, 1] — > R by 



for all 9 E (0, 1). From the assumptions we know h(X) = 0 and h(l) = 0. Since h is C 1 on [A, 1], by 
the mean- value theorem there must exist 9 E (A, 1) such that h'(9) = 0, which is a contradiction. 

For the purposes of this paper, we would like to construct a convex set of functions whose sub-level 
sets are connected. Although the convex set of all convex functions on M n will suffice, using it would 
unnecessarily restrict the class of sets describable to be convex. One cannot simply use the set of all 
functions whose 0-sub-level set is connected, since this set of functions is not convex. We therefore 
choose the set of strictly star-shaped functions, which is a convex set. We will use strictly star-shaped 
polynomials to represent sets. This is significantly more general than existing approaches using quadratic 
functions [1-3]. Also, it has been shown in Lemma 3 that if g is strictly star-shaped, then sub(g) is strictly 
star-shaped. By using this property, we can easily pose the star-shaped constraints on g to make sub(g) 
a connected set. 

3.5. An Algorithm for Backward Advection 

Here we will state the result of the backward advection algorithm. Given a strictly star- shaped 
polynomial gr i _ 1 such that sub^g^i) C R, and sub(gi-i) is bounded and positively invariant, we compute 
a polynomial ^ such that sub(Ahgi) ~ sub(gi-i) as follows. 

Pick a > 0 and positive integer d. Solve, using semidefinite programming, the following feasibility 
problem for g { E M[x], si, s 2 , s 3 , s 4 E E. 



h(9) = g(6y) 



for 9 E [0, 1] 



Then the derivative of h is 



h'{9) = l -Dg{9y){9y) 



> 0 



9i(0) = ~ 
Dgi(x)x > 0 

S3 - H9i-i + B {h _ a) gi = 0 
si + s 2 gi-i - B h gi = 0 
deg(ft) < d 



1 



Here we introduced an important parameter, a, which we think of as follows. The above algorithm 
finds a degree d polynomial g { such that g { is strictly star shaped, 4>y l sub(g i ) C sub(g i ^ l ), and 
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(fih-a sub(gi) D sub(gi-i). Hence the parameter a may be thought of as a tolerance that allows for the 
constraint that gi is required to have degree d or less. Then from the result of Theorem 1, lim^oo sub(gi) 
converges to the domain-of-attraction. It should be noted that this technique only works in the case that 
the advected set is positively/negatively invariant. 

3.6. Stopping Conditions 

By using the proposed level-set method, one can successfully propagate the system states backward in 
time. However, a stopping criterion is still needed to terminate the iterations. To detect the convergence 
of the advected sets, the closeness of two semi-algebraic sets is analyzed. The following result shows 
that the closeness of two semi-algebraic sets can be estimated by using scaled sets. The detailed proof 
can be found in [10]. 

Theorem 2. Suppose g\ and g 2 are strictly star- shaped functions, sub(gx) C sub(g 2 ), and sub(g 2 ) is 
bounded. Suppose x\, x 2 G M. n are two points such that 

x\ G d sub(gi) and x 2 G dsub(g 2 ), 

and X\ = ax 2 for some a > 0. Define the function q : R n — > R by 

q(x) = g 2 {\x) for all x G R n 

where A > 1. Then ifsub(q) C sub(gi), 

fe^ili<l_ A -i 

IrjII 

Figure 2. Example of stopping conditions. Curve 1 is the inner set and curve 2 is the outer 
set. Curve 3 is the shrunk set of the outer set. 

1.5 

1 

0.5 

0 

-0.5 
-1 

"-1.5 -1 -0.5 0 0.5 1 1.5 

To determine when the algorithm should terminate, one formulates an optimization problem using 
Lemma 1 to determine the smallest A > 1 such that sub(q) C sub(g 2 ), where q(x) = g 2 (Xx). Again, 
this may be evaluated using semidefinite programming. In practice, one picks a A > 1 in advance, 
and checks this condition after each iteration. Figure 2 shows an example. Here Curve 1 is dsub(gi), 
Curve 2 is d sub(g 2 ), and Curve 3 is d sub(q). The largest radial deviation between Curves 1 and 2 is 
less than 0.3. 
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3. 7. The Local Lyapunov Function 

We find a local Lyapunov function in order to construct an initial star-shaped positively invariant set. 
The following result is standard. 

Proposition 1. Suppose f : W 1 — > M. n is analytic and the origin is a stable equilibrium point. Also 
suppose V : W 1 — > K is a C 1 function, 7 > 0, and the set 

D 1 = {xeR n \ V(x) < 7} 

15 compact. Further suppose 

V(x) > 0 for allx=£0 
V(0) = 0 

DV(x)f(x) < 0 forallx^0,x E D y 
Let g 0 (x) = V(x) — 7. Then sub(g 0 ) is positively invariant, and sub(g 0 ) C R. 

One simple approach to finding an initial sub-level set is to find a quadratic Lyapunov function for the 
linear model of the system, and use a small sub-level set of this quadratic polynomial as the initial set. 

An alternative method which often gives a much larger initial set is as follows. Choose a polynomial 
p G M[x) such that sub(p) C R. We then solve the following convex feasibility problem. Find V G M[x] 
and s 0 ,si e £ such that 

DV(x)x>0 forallx^O 
V{x) > 0 for all x ^ 0 
V(0) = 0 

DV(x)f(x) + s 0 - sip = 0 for all x ^ 0 

Similar methods for finding local Lyapunov functions along with details on the construction of the 
associated semidefinite program may be found in [5,8]. Here we have added the first constraint to ensure 
that V — 7 is strictly star-shaped for 7 > 0. Note that these constraints imply that all sub-level sets of V 
are compact. Given V, we then solve the convex program 

maximize 7 

subject to V — 7 — s 0 — s x p — e = 0 for all x 

Sq, Si G £ 

where e > 0 is small. The optimal 7 satisfies sub(V — 7) C sub(p). Then V and 7 satisfy the assumptions 
of Proposition 1 and so we may use g 0 = V — 7 as the function defining our initial level-set. 

After we have found an estimate of the DoA, we can then use Proposition 1 to find a Lyapunov 
function that can be used to describe the behavior of the system within the DoA. To adequately describe 
the system behavior, this approach requires that we have a good estimate of the DoA. The following 
examples show that the level-set algorithm can precisely estimate the DoA. 
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3.8. Examples of Do A Estimation 

Example 1. Consider the following dynamical system 

x = 0.5y - x(l - x 2 - 0.25y 2 ) 
y = - x - o.5y(l - x 2 - 0.25y 2 ) 

The origin is a locally stable equilibrium point. Here we start with the initial polynomial 
go = 2x 2 + 2y 2 — 1. The results of the level-set method are shown in Figure 3, using time step h = 0.2. 
It can be seen that the successive iterates approach the true DoA. 

Figure 3. Successive iterates of the level-set algorithm. 

2 
1.5 
1 

0.5 
0 

-0.5 
-1 
-1.5 

~-2 -10 12 

Example 2. Consider the Van der Pol oscillator with inverted time: 

x = -y 

y = x-y(l- x 2 ) 

This system is locally stable around the origin. An initial sub-level set given by the quadratic 
polynomial g 0 = Ax 2 + Ay 2 — 1 is used which can be verified to be positively invariant. A time step 
ofh = 0.2 is used. The time tolerance parameter a is 0.02. The even-numbered iterates go, g^, #4, . . . 
are shown in Figure 4. For reference, some of the iterates are listed below and normalized to allow 
integer coefficients. 

p 2 = -1, 000 + 2, 252 y 2 - 88 y 4 + 11 y 6 - 907 xy - 56 xy 3 - 4 xy 5 + 3, 883 x 2 
+ 360 x 2 y 2 - 57 x 2 y 4 + 660 x 3 y - x 3 y 3 - 417 x 4 + 21 x 4 y 2 + 81x 5 y + 260 x 6 

Pi = -1,000 + 1, 6Uy 2 - 137 y 4 + 16y 6 - 1, QhAxy - 170 xy 3 + \Axy b + 3, 162 x 2 
+ 480 x 2 y 2 - 43 x 2 y 4 + 94 x 3 y - 35 x 3 y 3 + 144 x 4 - 2 x 4 y 2 + 192 x 5 y + 335 x 6 
p 2S = -10, 000 + 2510 y 2 - 56 y 4 + 2 y 6 - 4, 306 xy + 42 xy 3 + 4 xy 5 + 4, 099 x 2 
+ 25 x 2 y 2 + 2 x V + 1, 103 x 3 y - 27 x 3 y 3 - 687 x 4 - x 4 y 2 + 2x 5 y + 84 x 6 

It can be seen that the iterates gradually approach the exact boundary of the DoA. After 30 iterations, 
the solution covers most of the stable region. 
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After 40 iterations, the stopping criteria allowing an absolute radial change of 0.01 has been 
met. The final result is shown in Figure 4 as Curve 4. Curve 1 is dC(go) if we were using the 
semidefinite -programming based procedure in Section 3.7. For comparison, Curve 2 is the result of [2] 
and Curve 3 is the result of [1 ]. 

Figure 4. Van der Pol oscillator. The left figure shows the sequence of iterations. The right 
figure shows the final iterate as Curve 4 along with some other results. Curve 1 is d sub(g 0 ), 
when g 0 is obtained through the semidefinite programming based procedure in Section 3.7. 
For comparison, Curve 2 is the result of [2] and Curve 3 is the result of [1]. 




Example 3. The following dynamic system is the Example S4 taken from [3]. 

x = -2x + y + x 3 + y 5 
y = —x — y + x 2 y 3 

Here the initial set is a small circle around the origin. After 20 iterations, the estimated boundary of the 
DoA has reached the pre-specified bound. The final iterate is shown in Figure 5. The dashed curves show 
the two system trajectories which are used to represent the true boundary of the domain of attraction. 
Curve 2 represents the result of the final iterate. Curve 1 is the result from [3 ]. 

Figure 5. Sub-level sets of go and g 2 o- The dashed-dot curve is the result from [3] and the 
solid curve is the result of the advection algorithm. 
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4. Configuration of a PLL 

Figure 6 shows the basic configuration of a PLL. It has three components; a phase detector, a loop 
filter, and a voltage controlled oscillator(YCO). The VCO generates an output signal whose phase, 
9o(t), depends on the phase, of the input signal. The PLL is phase locked when the phase error 
4>(t) = 9i(t) — 9 0 (t) is a constant value and the loop is in stable equilibrium state. Usually, it is desired 
that the phase error, 4>{t), is maintained at zero. 

Figure 6. Basic Configuration of a PLL. 



e, 



Phase 


4> 


Low Pass 




Voltage 
Controlled 
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— ► 


Filter 


► 



e 



Of interest is the behavior of the phase error 4>{t). Because of its sinusoidal nonlinearity in the PLL, 
the phenomenon of chaos is believed to exist [11,12] and its inherent chaotic behavior for broadening the 
pull-in range of PLL has also been realized [13,14]. A nonlinear controller can drive PLL from chaotic 
state into periodic state or vice versa [15]. For higher-order PLL, it is not possible to determine whether 
the loop will or will not slip cycles using the initial frequency alone. In this case, one might define the 
pull-in range as the separatrix ordinate at 0 = 0 [17]. Analyzing the DoA of the PLL system provides 
a better description of the region in which a PLL locks up without slipping. The Lyapunov method has 
been used for stability analysis in control systems. Here the advection algorithm will be used to find 
the guaranteed stability boundary of the PLL system and the associated local Lyapunov function is then 
used to further refine the controller parameters. In [16], a Lyapunov styled analysis for PLL system up 
to third order is presented. The method shown in this section provides a way to analyze the DoA for a 
more general system. Also, the form of the Lyapunov function used here is much more flexible. 

Figure 7 shows the nonlinear model of the PLL. The sine function here represents the phase detector 
of the system. K in Figure 7 stands for the loop gain of the system. F(s) is equivalent to the low pass 
filter shown in Figure 6 and it corresponds to the controller of the PLL. Finally, the integrator in Figure 7 
is the voltage or numerically controlled oscillator. The key idea of a PLL system is to use the command, 
1/2, from F(s) to steer the oscillator such that 6o(t) tracks 0»(t) as closely and quickly as possible. 



Figure 7. Model of the Phase-Locked Loop. 
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4. 1. Second Order PLL 



To use the nonlinear design approach, start with a reference design. The reference design used in this 
paper is the linear model of a PLL system. A Proportional-Integrator (PI) controller is chosen to be the 
filter, F(s), as 

F(s) = i±™ (3) 

Using the model shown in Figure 7, it is routine to check that the resulting dynamic equation of the 
system is 

— + K— cos {(f)) — + —K sm 0 = — (4) 
dt 1 t\ dt T\ dt 2 

Assume that the received signal frequency is varying linearly with time and has zero radial acceleration 

and let x x = 0, x 2 = 4>. The PLL system can be rewritten as the following state space model 

X\ = X2 

T 2 1 

x 2 = —K — cos(xi)x 2 i^sin(a;i) (5) 

= ki cos(xi)x 2 + k 2 sin^x) 

Equation (5) is the final nonlinear model of the second order PLL. A linearized model can then be 
derived as 



X\ = X2 

x 2 = hx 2 + k 2 xi 

The filter F can then be designed using existing linear design approach [17]. One typical choice is to 
let uj n = 15, C = 0.707, K — 1, where u n is the natural frequency, £ is the damping ratio, and K is the 
overall gain. The two coefficients of the filter are then given as T\ = r 2 = ^f- 

5. Phase-Locked Loop Analysis and Design 

In this section, the second order PLL controller will be used to demonstrate the nonlinear design 
approach. The same approach can also be applied to the third order controller design. 

5. 1. Pull-In Range of the Traditional PLL System 

Now the advection algorithm can be applied to the PLL nonlinear system. To reduce the required 
number of iterations, a local Lyapunov function is used as the initial set. After a few iterations of 
the algorithm, it gives us the estimated domain-of-attraction of the system. The result is shown in 
Figure 8. Note that the estimated region is based on the Taylor series expansion of the sine and cosine 
functions. From Ston-Weierstrass Theorem, we can approximate the sine and cosine functions to the 
desired accuracy within a bounded interval. In this example, two degree- 10 polynomials are used to 
approximate sine and cosine functions and the estimated region is only valid between — n and tt. More 
terms of the Taylor series could also be used to improve the accuracy. 



Sensors 2011, 11 



6588 



Figure 8. Result of the original PLL system. 
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5.2. PLL System Controller Design 

After getting a good estimated DoA, a local Lyapunov function that describes the system behavior 
can be easily computed using Proposition 1 . This local Lyapunov function can also be used to describe 
the trajectories of the system in the DoA. Figure 9 shows several Lyapunov level-sets obtained from the 
SOS approach. 

Figure 9. Local Lyapunov level sets of the original PLL system. 




-4-3-2-101234 



The SOS techniques can be applied to the design of a controller. For the PLL system, it is desired 
to design a system which has a larger DoA or faster converging speed. Here, the objective is to find 
possible system parameters such that the same Lyapunov function is still valid and the system converges 
faster. This can be done by solving a semidefinite program. 
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Suppose V is a local Lyapunov function for the PLL system. Using the SOS technique, solve the 
following optimization problem: 



where a, b G M+ specify the domain of the constraints and q is a positive definite performance 
polynomial specified by the user. p(k) is a linear constraint of controller parameter k. As before, 
Si, s 2 , s 3 , s 4 are SOS polynomials. 

Since V is now a given function, the above constraints are linear in the controller parameters. The first 
constraint shows that V is a valid Lyapunov function in sub(V — a). This constraint is used to specify 
the desired DoA to maintain. The second constraint along with the objective function will put an upper 
bound on the derivative of the Lyapunov function in sub(V — b). Faster decreasing rate implies faster 
converging speed. This specifies the performance requirement of our system. The user could also put 
different performance requirements in different sub-level sets of V. 

Besides dynamic performance constraints, noise bandwidth constraints will be applied as well. The 
noise bandwidth for this PLL system with PI controller has the following form [17] 



Assume £ > 1, u n > 1. Then 



This linear upper bound will be used to find a set of controller parameters that have better dynamic 
performance while maintaining the same noise bandwidth. 



max a 



s.t. - (DV)f = si + s 2 (a - V) 

~(DV)f-aq = s 3 + s,(b-V) 
p{k) < 0 



(6) 




Figure 10. Comparison of the domain-of-attraction. Left: linear design. Right: nonlinear 
design. 
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The system phase portrait as well as the estimated stable region are shown in Figure 10. The nonlinear 
design has K = \, u n = 10.813, and £ = 1.3303. The noise bandwidth is 8.1082 Hz, which is slightly 
higher than the noise bandwidth of the linear design, 7.9546 Hz. From the phase portrait, the nonlinear 
design has approximately 20% larger guaranteed domain-of-attraction. It can also be observed from 
the phase portrait that the nonlinear design has less overshoot than the linear design. This shows that 
the proposed method increases the system performance while not sacrificing too much of the noise 
rejection capability. 

A Simulink model is used to compare the nonlinear designed PLL controller with the linear design. In 
this Simulink model, the sinusoidal input is collapsed by measurement noise and clock noise with zero 
mean and variances 0.1 and 0.0001, respectively. Figure 11 is the phase error of the two designs. It is 
clear that the nonlinear designed controller has a much faster convergence rate than the original design. 

Figure 11. Phase error of Simulink simulation. Left: linear design. Right: nonlinear design. 





Figure 12. Real GPS experimental results. Left: linear design. Right: nonlinear design. 
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Both systems are tested on a NORDNAV R25 software GPS receiver. This receiver collects, 
down-converts and samples the GPS data by the front end, so that the collected GPS data can be 
post-processed repeatedly using different tracking-loop filter orders. The results are shown in Figure 12. 

It can be seen that the original system has some overshoot and converges around 300 ms. The 
nonlinear design has much less overshoot and converges about five times faster than the original design. 

6. Summary 

In this paper, we presented a method of designing a PI controller of a PLL system. This design 
approach is based on the polynomial nonlinear model of the PLL system. This approach starts with 
the linear design of the controller and then estimates the DoA of the linear designed system to get the 
suitable local Lyapunov function for the system. The Lyapunov function is then used as the performance 
constraints to further refine the performance of the system outside the linear region. The DoA of the 
initial design can also be extended to get a better pull-in region of the PLL system. This approach gives 
us a way to design a fixed form controller for a nonlinear system. 
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